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Abstract 

We investigate the conditions for reactivity enhancement of catalytic processes in porous solids by 
use of molecular traffic control (MTC) as a function of reaction rate and grain size. With dynamic 
Monte-Carlo simulations and continuous-time random walk theory applied to the low concentration 
regime we obtain a quantitative description of the MTC effect for a network of intersecting single-file 
channels in a wide range of grain parameters and for optimal external operating conditions. The 
efficiency ratio (compared with a topologically and structurally similar reference system without 
MTC) is inversely proportional to the grain diameter. However, for small grains MTC leads to a 
reactivity enhancement of up to approximately 30% of the catalytic conversion A — > B even for 
short intersecting channels. This suggests that MTC may significantly enhance the efficiency of a 
catalytic process for small porous nanoparticles with a suitably chosen binary channel topology. 
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I. INTRODUCTION 



Zeolites are used for catalytic processes in a variety of applications, e.g. cracking of large 
hydrocarbon molecules. In a number of zeolites diffusive transport occurs alongquasi-one- 
dimensional channels which do not allow guest molecules to pass each other Due to 
mutual blockage of reactand A and product molecules B under such single-file conditions 
P| the effective reactivity of a catalytic process A — ► B - determined by the residence time 
of molecules in the zeolite - may be considerably reduced as compared to the reactivity in 
the absence of single-file behaviour. It has been suggested that the single-file effect may be 
circumvented by the so far controversial concept of molecular traffic control (MTC) (3, |j|- 
This notion rests on the assumption that reactands and product molecules resp. may prefer 
spatially separated diffusion pathways and thus avoid mutual suppression of self-diffusion 
inside the grain channels. 

The necessary (but not sufficient) requirement for the MTC effect, a channel selectivity 
of two different species of molecules, has been verified by means of molecular dynamic (MD) 
simulations of two-component mixtures in the zeolite ZSM-5 \^ and relaxation simulations 
of a mixture of differently sized molecules (Xe and SF 6 ) in a bimodal structure possessing 
dual-sized pores (Boggsite with 10-ring and 12-ring pores) (q). Also equilibrium Monte-Carlo 
simulations demonstrate that the residence probability in different areas of the intracrys- 
talline pore space may be notably different for the two components of a binary mixture [3] 
and thus provide further support for the notion of channel selectivity in suitable bimodal 
channel topologies. 

Whether a MTC effect leading to reactivity enhancement actually takes place was ad- 
dressed by a series of dynamic Monte Carlo simulations (DMCS) of a stochastic model 
system with a network of perpendicular sets of bimodal intersecting channels and with cat- 
alytic sites located at the intersecting pores (NBK model) P, llO, lllfl . The authors of these 
studies found numerically the occurrence of the MTC effect by comparing the outflow of 
reaction products in the MTC system with the outflow from a reference system with equal 
internal and external system parameters, but no channel selectivity (Fig. [TJ. The depen- 
dency of the MTC effect as a function of the system size has been investigated in [l^j]. The 
MTC effect is favored by a small number of channels and occurs only for long channels be- 
tween intersections, which by themselves lead to a very low absolute outflow compared to a 
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FIG. 1: REF system (left) with N = 3 channels and MTC system (right) of the the same size. In 
contrast to the REF case, where we allow both types of particles (A and B particles) to enter any 
channel, in the MTC system A particles are carried through the vertical a channels whereas the 
B particles diffuse along the horizontal (3 channels. Black squares indicate catalytic sites where a 
catalytic transformation A — > B is allowed. 



similar system with shorter channels. A recent analytical treatment of the master equation 
or this stochastic many-particle model revealed the origin of this effect at high reactivities 
~\\ . It results from an interplay of the long residence time of guest molecules under single- 
file conditions with a saturation effect that leads to a depletion of the bulk of the crystallite. 
Thus the MTC effect is firmly established, but the question of its relevance for applications 
remains open. 

Here we address this question by a systematic study of the MTC effect as a function of 
the reactivity of the catalytic sites and as a function of grain size, but using fixed small 
channel length. This choice is motivated by potential relevance from an applied perspective. 
Moreover, for the same reason we determine the MTC effect by making a comparison with 
the reference system using the same set of fixed internal (material-dependent) parameters, 
but (unlike in previous studies [ill 12, 13) for each case (MTC and REF resp.) different 
optimal external (operation-dependent) parameters which one would try to implement in 
an industrially relevant process. It will transpire that a significant MTC effect (reactivity 
enhancement up to ~ 30%) occurs in our model system even for small channel length at 
realistic intermediate reaction rates of the catalyst, provided the grain size is sufficiently 
small. This may be of interest as since the first successful synthesis of mesoporous MCM-41 
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nanoparticles [IJ], there has been intense research activity in the design and synthesis of 
structured mesoporous solids with a controlled pore size. In particular, synthesis of bimodal 
nanostructures with independently controlled small and large mesopore sizes has become 
feasible 

In this work we do not study a specific process in a specific material, but we demonstrate the 
validity of the MTC concept even if channels inside the porous material are short. This is 
novel and - from an applied viewpoint - crucial since it is a necessary condition for starting 
expensive and time-consuming quantitative investigations in specific settings. 



II. NBK MODEL 

As in [12I, [3] we consider the NBK lattice model j$] with a quadratic array of N x N 
channels (Fig. [TJ which is a measure of the grain size of the crystallite. Each channel has L 
sites between the intersection points where the irreversible catalytic process A — > B takes 
place. We assume the boundary channels of the grain to be connected to the surrounding gas 
phase, modelled by reservoirs of constant densities such that the entrances of the respective 
channels (extra reservoir sites) have a fixed A particle density p. We assume the reaction 
products B which leave the crystallite to be removed immediately from the gas phase such 
that the density of B particles in the reservoir is always 0. Short-range interaction between 
particles inside the narrow pores is described by an idealized hard core repulsion which 
forbids double-occupancy of lattice sites. 

The underlying dynamics are stochastic. We work in a continuous time description where 
the transition probabilities become transition rates and no multiple transitions occur at the 
same infinitesimal time unit. Each elementary transition between microscopic configurations 
of the system takes place randomly with an exponential waiting-time distribution. Diffu- 
sion is modelled by jump processes between neighbouring lattice sites. D is the elementary 
(attempt) rate of hopping and is assumed to be the same for both species A, B of particles. 
In the absence of other particles D is the self-diffusion coefficient for the mean-square dis- 
placement along a channel. If a neighboring site is occupied by a particle then a hopping 
attempt is rejected (single-file effect). The dynamics inside a channel are thus given by the 
symmetric exclusion process [3, llTl . fisl which is well-studied in the probabilistic [^J 
and statistical mechanics literature |2lj. The self-diffusion along a channel is anomalous, 
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the effective diffusion rate between intersection points decays asymptotically as 1/L, see [3] 
and references therein. 

At the intersections the reaction A — > B occurs with a reaction rate c. This reaction rate 
influences, but is distinct from, the effective grain reactivity which is largely determined 
by the residence time of guest molecules inside the grain which under single-file conditions 
grows in the reference system with the third power of the channel length L 2^|. At the 
boundary sites particles jump into the reservoir with a rate D(l — pa — Pb) in the general 
case. Correspondingly particles are injected into the grain with rates DpA,B respectively. 
As discussed above here we consider only pa = p, Pb = 0. 

For the REF system A and B particles are allowed to enter and leave both types of 
channels, the vertical (a) and horizontal ((3) ones. In case of MTC A(B) particles will 
enter a(/3)-channels only, mimicking complete channel selectivity. Therefore all channel 
segments carry only one type of particles in the MTC case. For the boundary channels 
complete selectivity implies that a-channels are effectively described by connection with an 
^4-reservoir of density pa = p (-B-particles do not block the boundary sites of a-channels) 
and /3-channels are effectively described by connection with a 5-reservoir of density ps = 0, 
respectively. (A-particles do not block the boundary sites of /3-channels.) This stochastic 
dynamics, which is a Markov process, fully defines the NBK model. 

In both cases, MTC and REF system, the external concentration gradient between A and 
B reservoir densities induces a particle current inside the grain which drives the system in a 
stationary nonequilibrium state. For this reason there is no Gibbs measure and equilibrium 
Monte-Carlo algorithms cannot be applied for determining steady state properties. Instead 
we use dynamic Monte-Carlo simulation (DMCS) with random sequential update. This 
ensures that the simulation algorithm yields the correct stationary distribution of the model. 

III. MONTE CARLO RESULTS 

Anticipating concentration gradients between intersection points we expect due to the 
exclusion dynamics linear density profiles within the channel segments 12, 2lj], the slope 



and hence the current being inversely proportional to the number of lattice sites L. The 
total output current j of B particles, defined as the number of -B-particles leaving the grain 
per time unit in the stationary state, is the main quantity of interest. It determines the 
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FIG. 2: juEF (solid symbols) and jure (open symbols) as a function of the reservoir density for a 
system with iV = 5, L = 1, c = 0.1. 

effective reactivity of the grain. 

We are particularly interested in studying the system in its maximal current state for 
given reactivity c and size constants N, L, which are intrinsic material properties of a given 
grain. The A particle reservoir density p, determined by the density in the gas phase, can be 
tuned in a possible experimental setting. Let us therefore denote the reservoir density which 
maximizes the output current by p* and the maximal current by j*. For MTC systems as 
defined above we always expect p* MTC = 1, since the highly charged entrances of a-channels 
do not block the exit of -B-particles and hence do not prevent them from leaving the system. 
Fig. 01 shows j as a function of p for both a MTC and REF system of N = 5, L = 1 and 
reactivity c = 0.1. Indeed for MTC the maximal output occurs for the maximal reservoir 
density. In case of the REF system p* REF as well as j* REF need to be found by simulation. 
We iteratively approach the maximal current by a set of 9 datapoints. The "best" datapoint 
has been chosen in order to approximate the maximum. Statistical errors are displayed. 
They are, however, mostly within symbol size. 

In order to measure the efficiency of a MTC system over the associated REF system we 
define the efficiency ratio 

Jmtc 



R(c, N, L) 



Jref 



(1) 



which is a function of the system size N, L and reactivity c. 

Fig. El (left) shows the measured ratio R for a large range of reactivities. We plotted 
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j MTC /j REF for L=2 and N=3..1 5 - KBN geometry 



Maximal Ratio j MTC /j REF L=1,2,3 (KBN) 



1.8 
1.6 
1.4 
1.2 
cc 1 
0.8 
0.6 
0.4 
0.2 







m 

o 






- 


i I 




I 










«.... 


m : 
m 










■ e 




N=3 






















■ 

o ■ ■ 


N 5 










* I • • . 












;. j «. 

7 .t .* S . 


. N=15 





0.0001 0.001 0.01 0.1 1 10 100 1000 



1.2 

1 

0.8 



















i 


S \ 












- 




L=3o 

















■ i 

L=2 
: ■ 








o 






m 


5 L-1 




■ 




1 








V m 


□ 


□ 


i 






□ 




i 















2 4 6 



10 12 14 16 18 20 
N 



FIG. 3: Left: Ratio R(c) for different number of channels N and L = 2. Right: Maximal ratio R* 
for different L 

systems with L = 2 and different N. We note that the MTC effect has a strong negative 
dependence on increasing N for all N and increasing c above some optimal value c*. We 
denote this optimal value by R ma x- Fig. [HI (right) shows R* max {N) for different L and proves 
that there is an MTC effect for any L and any N. Notice, however, that with increasing N 
the optimal ratio not only decreases, but appears at unnaturally small reactivities c. This 
is highly undesirable as then the actual output current shrinks to zero. Even though more 
efficient than a REF system, a MTC grain with large N would not operate under practically 
relevant conditions. 

As expected from previous results Q the MTC effect is seen to increase with in- 
creasing L, even though the measure R used here is different (Fig. HJ). This follows from 
theoretical studies of single-file systems. The mean traveling time of a product molecule 
through a channel of length L is proportional to L 2 in the MTC case as in ordinary diffu- 
sion, but proportional to L 3 in the REF case due to mutual blockage. Hence the current 
is proportional to 1/L in a MTC system, but proportional to 1/L 2 in a REF system. This 
holds for all values of the parameters and hence for sufficiently large L the MTC system 
becomes more efficient. 

In order to understand the strong decrease of R for large reactivities and large number 
of channels it is instructive to study the stationary density profiles. Fig. El shows the B 
particle densities for a REF system (left) and MTC system (right) with N = L = 9 and 
large reactivity c — * oo. Theoretical investigations of MTC systems in the large reactivity 
case (3| show that in the state of maximal current (p = 1) the output of B particles is 
independent of the number of channels. For large L the maximal output current becomes 
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j MTC /j REF for L=1 ,2,3 and N=1 5 (KBN) 
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FIG. 4: Dependence of the output ratio on the channel length L. 

J*->oo — wnere D is the diffusion constant. A nonvanishing current of B particles can be 
observed only at the four corners of the lattice (Fig. El right). 

For fixed moderate c this extreme situation is not realized, but nevertheless with increas- 
ing N one expects that the bulk gets increasingly depleted, since in each layer a fraction 
of A particles is converted into B particles. Thus the total A-density in each layer may be 
described by the form 



ax 



-jN A (x) 



(2) 



The coarse-grained ansatz for the number N A (x) of A-particles in layer x predicts an expo- 
nential decrease of the A density, leaving only an active boundary layer of finite thickness 



£ = I/7 oc 1/c 



(3) 



at the top and bottom respectively of the (in our simulation two-dimensional) grain. Hence, 
as a function of N, j* MTC saturates at some constant 



lim j* MTC (c, N, L) = C* MTC {c, L). 

N —>oc 



(4) 



On the other hand, in the REF system the output current scales linearly with increasing 
N for all, even large, c (Fig. El left). This is because even though again the bulk depletes 
with increasing N the active boundary layer is a surface scaling linearly with N. Thus 



Jim Jref(c, N, L) = NC* REF (c, L) 



(5) 
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Profile of B particles - REF {N=L=9, c=lnf) Profile of B particles - MTC (N=L=9, c=lnf) 




FIG. 5: Profiles of the REF (left) and MTC (right) system in the large-reactivity case. N = L = 9 
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FIG. 6: Active channel segments in the large-reactivity case. REF system (left), MTC system 
(right) 

Hence 

R(c, N, L) oc l/N (6) 

and the MTC effect vanishes at some N for fixed reactivity c and channel length L. 

Looking on a wider range of reactivities (Fig. EJ) we notice that the maximal currents for 
MTC systems saturate for some intermediate c whereas in the REF case the current reaches 
a plateau for large reactivities. This observation can be rationalized by noticing that an 
increase of the output with increasing c is limited by the incoming current of available A 
particles. Since in the REF system A and B particles block each other an increasing current 
of -B-particles always restricts the number of A particles diffusing in. Hence the saturation 
due to high reactivity sets in only for larger values of c than in the MTC system. 
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FIG. 7: Maximal currents for the REF system (left) and the MTC system (right). L = 2 and 
different N. 

IV. REF WITH SMALL REACTIVITIES 



The arguments put forward above for explaining qualitative and quantitative features of 
the MTC effect at intermediate and large reactivities c fail for small reactivities, i.e., when 
c is of the order of the inverse mean intracrystalline residence time. We first consider the 
reference system. 

In this case the grain contains only a very small number of B particles at any time. In this 
low-concentration regime the diffusion of B particles inside the grain may be described by a 
linear diffusion equation with an effective diffusion coefficient determined by the interaction 
with A particles as follows. The A particles are considered as a medium with constant 
equilibrium density p throughout the system. At the intersections B particles are created 
randomly with effective rate cp. We then describe the diffusion of a B particle in the channel 
between intersections by a random walk from intersection to intersection with an effective 
diffusion constant -D e // given by the (inverse) mean travelling time between intersections. 
Let P( x ,y) be the B particle density at the intersection denoted by (x, y) and A a discrete 
two-dimensional Laplace operator 

1 



^P(x,y) = ^ (p(x+l,y) + P{x-l,y) + P(x,y+l) + P(x,y-1) ~ ^P(x,y)) ■ 



(7) 



Here the lattice unit is given by the channel length rather than the pore size inside a channel. 
For the B particle density at intersections we thus obtain 

d 

Qj.P{*,V) = D eff^P(x,y) + Cp. (8) 

A stationary solution to (JHJ) can be found by use of the discrete sine transform pt q> p) 
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J2x=i ~52y=i P(*,y) s * n ~nTi s * n N+T- We ex P ress © (with vanishing time derivative) in terms 
of the transformed density P( q , p ) ■ Taking into account the boundary conditions po iV = p x $ = 
Px,n+i = Pn+i, v = with 0<x, y < N + 1 we find 

2cp A v^v^ ■ Q nn PKm , 

««■> = C037fe + COS7gT - 2 E E SU1 jvTT sm nTT- (9) 

iV+l iV+1 n= l JTl=l 

The non zero contributions of the double sum can be expressed as a product of two Cotan- 
gents. Transforming back finally yields 

N N D 

" M = " £ £ cos ^ + cos ft - 2 Sm NT-l Sm WTT 




if n or m even 

-B(n,m) = { (11) 



_ rot m?r rot n7T else 
(AT+1) 2 L,UL 2(N+1) L ' UL 2(jV+l) elSt; 



-D e // together with the reactivity and reservoir density, is absorbed into a fitting parame- 
ter A ~ jjjj- With pfx™) being the B particle densities obtained from simulations and pf^y) 
the theoretical densities, we define the homogenized mean square deviation 



n Sim rt^^ 1 

E [ XIS? ] w 



\ intersections V (z,2/) H {x,y) , 



as a measure for the quality of the approximation. The sum runs over all intersections, 
with the local deviation normalized by the local mean. This assures that all intersections 
contribute with their proper weight. 

For large c the profile flattens (see Fig. EJ) and is not very well described by (JT0| . Fig. [HI 
shows Q as a function of c. The boundary density has been chosen to p = 0.5. For small 
reactivies the collapse of the simulated and calculated profile is fairly good (Q ~ 0.1). 
However, the best collapse occurs for small intermediate reactivities. This somewhat 
peculiar behavior can be explained by assuming that this random walk is not a Markov 
process as implied in the derivation of (fT0|) . The structural change in the occupation of the 
channel segment, once a particle covered the distance between two adjacent intersections, 
implies a "memory" effect. Hence we cannot assume perfect statistical independence of 
subsequent random traveling times between intersections which implies that the Markov 
assumption is not very well satisfied. As we increase the reactivity more than one B 
particle may be present in the system. This corresponds to an ensemble average over 
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FIG. 8: Collapse of the simulated and calculated profiles for REF systems. Left: N = 5, p = 0.5 
and different L. (right) N = 5, L = 1 and different p. 

almost independent random walkers which cancels the time correlations of each individual 
walker and thus improves the validity of the Markov assumption. This leads to the good 
data collapse for some intermediate reactivities. As c increases further both the assumption 
of equilibrium of A particles and the low-concentration approximation (JHJ) for B particles fail. 



V. MTC WITH SMALL REACTIVITIES 

Fig. El shows the stationary density profiles of a MTC system with small reactivity. We 
first note that the output current is proportional to the number N of (3 exit channels, in 
agreement with the observation that R ma x approaches a constant for large N, rather than 
decaying proportional to 1/N as R for fixed reactivity does. Moreover, due to rare transition 
events all a-channels are almost in equilibrium with the reservoir (Fig. [HI left). Thus it is 
sufficient to single out only one (3 channel. 

We adapt the approach we took for the REF system to the present case. Due to low 
concentration we use a linear diffusion equation for the 5-density inside a /3-channel. There 
are, however, two essential differences. (1) Because auf the absence of A-particles we do 
not consider the occupation of intersections alone, but we track the motion of B particles 
inside a channel. (2) Due to the inequivalence of the sites inside the /3-channels and the 
intersection points where also A may be located, we need to describe a random walker with 
space dependent hopping rates. Intersections, on the one hand, serve as sites of B particle 
creation with a rate proportional to the reactivity. On the other hand intersections are 
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Profile of A particles - MTC (N=L=5, c=0.0005) 



Profile of B particles - MTC (N=L=5, c=0.0005) 
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FIG. 9: Profiles for a MTC system N = L = 5 and small reactivity c = 0.005. (left) densities of A 
and (right) of B particles, p = 0.3 
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FIG. 10: (3 channel with hopping rates. L = 4 



occupied by A particles with a probability p and hence, block B particles. This leads to 
different hopping rates onto and from an intersection as displayed in Fig. EH 

The master equation description leads to a set of equations for the local B particle 
densities (n x ). 



d 



D (K-i> + (n x+1 ) - (n x ) - (1 - p) K» x = {L + l)r ± 1 
D(l-p)((n x ^) + (n x+1 ))-2D{n x )+cp x = (L + l)r (13) 
D ((n^i) + (n x+ i) - 2 (n,)) else 



Here x is the lattice position inside a channel in units of the pore size. In the stationary 
state left hand site of (fl"3|) vanishes and it follows a recurrence relation for the B-particle 
intersection densities p r = (^(L+i)r)- 



1 

Pr = 2 (Pr+i + + K ) 
K=^[{L + \)-(L-\)p\ 



(14) 
(15) 
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The inhomogeneity K depends on the segment length L and the transition rates. A 
solution satisfying (fT4~|) and the boundary condition po = Pn+i = can be obtained by use 
of a quadratic ansatz. We find 

p r = {n (L+l)r )= l -rK[N + l-r] (16) 

Due to the modified hopping rates in the neighborhood of the intersections the density 
profile is linear between (n^+iy-i) and (n(L+i) r -L) rather than between the intersections 
p r and pr-i itself. This becomes noticeable for large reservoir densities p. There is a "jump" 
between intersections and their adjacent lattice sites. Fig. [TT] illustrates this by means of 
a profile obtained from MC simulations). Let us consider a channel segment embedded 
by the two intersections p r and p r -i- Solving (fT3|) for the lattice sites located next to the 
intersections leads to 

/ \ Pr-l +pr(L-(L- l)p) 

- (L + l)-2Lp+(L-l)p2 (17) 

/ y Pr + Pr-l (L-(L- l)p) 

^ + D-l) ~ (L + l)-2Lp+(L-l)(? (18) 
and for the very first channel segment we find 

<•*> = rrrh/ < 19 > 

We have now fully determined the profile for MTC systems with small reactivity. The 
currents between neighboring lattice sites are proportional to the density difference ~ 
((n x +i) — ( n x)) with a proportional constant being the actual hopping rate. 

Fig. fT2l shows the simulated profile of a (3 channel together with the theoretical densities. 
We chose an intermediate reservoir density p = 0.5. The collapse of the two curves is very 
good. In order to study the quality of the approximation for different sets of parameter 
we use the definition (fT21) above. The normalized mean square deviation Q is plotted for 
different L and different p (Fig. [TJJ). In the regime of interest, i.e. for small reactivities, the 
profile is well described by the theory even for small L and rather high reservoir densities. 



VI. CONCLUSION 



Our simulations and analytical results describe the MTC effect quantitatively over a 
wide range of parameters within the NBK model. Moreover, trends which are independent 
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Profile of B particles - MTC (N=L=5, c=0.0005) 




FIG. 11: Demonstration of "jumps" at intersections for large p. Here we choose p = 0.9. 
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FIG. 12: Densities of a /? channel. Small open circles connected with the dotted line show the 
theoretical densities. Solid circles are densities obtained by MC simulation (N = L = 5, c = 0.0005, 
p = 0.5) . 

of model details have been identified by analytical calculations and lead to a more positive 
result than concluded by jl^] where no MTC effect at all was reported for short intercon- 
necting channels L — 1. For reasonable reactivities and channel lengths the MTC effect 
vanishes proportionally to 1/N, i.e., is inversely proportional to the grain diameter. This 
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FIG. 13: Collapse of the simulated and calculated profiles for MTC systems. Left: N = 5, p = 0.5 
and different L. Right: iV = 5, L = 5 and different p. 

was shown explicitly for a two-dimensional simulation model, but the reasoning that led to 
this conclusion extends straightforwardly to a three-dimensional system. Nevertheless, for 
optimized external process parameters the NBK model exhibits an enhancement of the ef- 
fective reactivity of up to approx. 30% for small grains and any (even short) channel length 
and reactivity c. In the present study we have not taken into consideration that smaller 
molecules may diffuse into larger channels, not only the smaller channels. This was done in 
order to keep the model simple so that the origin of the observed MTC effect becomes as 
transparent as possible. We did make simulations for a different model where this mecha- 
nism is taken into account. We have found in these simulations that the MTC effect becomes 
stronger. In summary, our investigations suggest that MTC may enhance significantly the 
effective reactivity in zeolitic nanoparticles with suitable binary channel systems and thus 
may be of practical relevance in applications. 
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